mfr_sensory = analyze_mfr_fmri(810001, 905000,.001);
mfr_imagery = analyze_mfr_fmri(906001, 1225000,.001);

figure(1);
plot(mfr_sensory,mfr_imagery,'.')
axis([0 140 0 140])
axis square

%figure(2)

N_neurons = length(mfr_imagery)

active_imagery = mfr_imagery>.5;
active_sensory = mfr_sensory>.5;

disp('% of all neurons active in imagery');
sum(active_imagery)/N_neurons*100
disp('% of all neurons active in sensory');
sum(active_sensory)/N_neurons*100

active = active_sensory|active_imagery;
N_active = sum(active)

disp('% of ACTIVE neurons active in imagery');
sum(active_imagery)/N_active*100
disp('% of ACTIVE neurons active in sensory');
sum(active_sensory)/N_active*100


disp('% of ACTIVE active only in imagery')
sum(active_imagery&~active_sensory)/N_active*100

disp('% of ACTIVE active only in perception')
sum(active_sensory&~active_imagery)/N_active*100

disp('% overlap');
overlap = active_sensory&active_imagery; 
sum( overlap )/N_active*100

disp('corr of imagery and perception population firing rates; ACTIVE neurons')
corr(mfr_imagery(active),mfr_sensory(active))

disp('corr of imagery and perception population firing rates; overlapping neurons')
corr(mfr_imagery(overlap), mfr_sensory(overlap))

figure(2)
boxplot([mfr_imagery(overlap), mfr_sensory(overlap) ]);


disp('ranksum mfr_imagery(active) vs mfr_sensory(active):');
[p h]=ranksum(mfr_imagery(active),mfr_sensory(active))
disp('median mfr sensory:');
median(mfr_sensory(overlap))

disp('median mfr imagery:');
median(mfr_imagery(overlap))


